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We investigate tiie diffusion coeiEcient of the time integral of the Kuramoto order parameter 
in globally coupled nonidentical phase oscillators. This coefficient represents the deviation of the 
time integral of the order parameter from its mean value on the sample average. In other words, 
this coefficient characterizes long-term fluctuations of the order parameter. For a system of 
coupled oscillators, we introduce a statistical quantity D, which denotes the product of N and the 
diffusion coefficient. We study the scaling law of D with respect to the system size A'^. In other 
well-known models such as the Ising model, the scaling property of _D is _D ~ Oil) for both coherent 
and incoherent regimes except for the transition point. In contrast, in the globally coupled phase 
oscillators, the scaling law of D is different for the coherent and incoherent regimes: D ~ 0(1/N°') 
with a certain constant a > in the coherent regime, and D ~ 0(1) in the incoherent regime. We 
demonstrate that these scaling laws hold for several representative coupfing schemes. 
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In real-world systems, large populations of coupled oscillators often experience global synchronous 
oscillations. To elucidate the general properties of such phenomena, considerable research has been 
conducted on simple models of globally coupled phase oscillators. The Kuramoto order parameter has 
been widely used to measure the degree of synchronization and to characterize the synchronization 
transition in the phase oscillator model. When the system size is infinite, fluctuations of the order 
parameter vanish after a transient period; the scaling law for this parameter has been well studied for 
a general coupling function. However, when the system is large but of finite size, the fluctuations in 
the order parameter do not vanish and their behavior has not been fully understood. It is not clear 
whether the conventional standard statistical quantities such as the variance and the correlation time 
of the order parameter can fully characterize its fluctuation behavior. Further, the dependence of the 
statistical properties of the order parameter on the coupling scheme is still not completely understood. 
As a step toward understanding these problems, we focus on a statistical quantity that characterizes 
long-term fluctuations in the order parameter. In other well-known models such as the Ising model, 
the scaling property of the statistical quantity with respect to the system size is the same for coherent 
and incoherent regimes except for the transition point. In contrast, in the globally coupled phase 
oscillators, the decay speed of the statistical quantity in the coherent regime is faster than that in 
the incoherent regime. This diflterence is caused by a diflterence in the correlations among the phases 
of the oscillators at diflterent times. We show that the scaling laws hold for a large class of general 
coupling schemes. 



I. INTRODUCTION 



Nonlinear systems are often used for modeling chemical reactions, engineering circuits, and biological populations 
Synchronization in such systems has attracted considerable attention in the past several decades. The phase 
description of the systems is one of the most effective methods to understand synchronization in interacting oscillatory 
systems P, 0] ■ Accordingly, there have been a number of studies on the globally coupled phase oscillator model Q , 
which is described as follows: 



K ^ 

(j = i,...,iv), (1) 

fc=i 



where Oj represents the phase of the jth oscillator, ojj represents the natural frequency of the jth oscillator, K > Q 
represents the coupling strength, h is the coupling function, and N is the number of oscillators. The oscillators are 
synchronized when the coupling strength is sufficiently large for the Kuramoto model where h{x) = sin(a;) |3|. The 
synchronization transition in the phase oscillator model ([IJ with infinite dimension (i.e., in the thermodynamic limit 
N — > oo) has been well studied with regard to its analogy to the second-order phase transition In particular, 

one of the main areas of focus in these studies has been the behavior of the order parameter R{t) as a measure of 
synchrony [21], which is defined as follows: 



Rit) = — 



N 



exp(27ri6'j) 



(2) 



A synchronization transition can be characterized by a change in the order parameter from zero to a non-zero value 
with an increase of K . We assume that the stationary state with R{t) = in the incoherent (desynchronized) regime 
supercritically bifurcates at the critical coupling strength K = Kc, above which the oscillators are synchronized 
or coherent. The scaling property of the order parameter around the synchronization transition point has been 
intensively studied: first, Kuramoto analytically investigated the phase oscillator model ([l} with the sinusoidal 
coupling function; then, Daido Crawford and Davies Q, and Chiba and Nishikawa [8,] considered the phase 
oscillator model ([1]) with more general coupling functions. 

However, finite size effects in the phase oscillator model ^ have not yet been fully understood. It is not clear 
whether the conventional standard statistical quantities such as the variance and the correlation time of the order 
parameter can fully reveal the characteristics of its fluctuations. Further, the dependence of the statistical properties 
of the order parameter on the coupling scheme is still not completely understood, because most previous studies on 
the finite size effects have only examined the Kuramoto model [9l-fl3| . It _is known that the scaling law of the order 
parameter in the infinite-size system depends on the coupling function [5|, |6|. However, it is not evident how the 
coupling function influences the scaling law of the order parameter in a finite-size system. Although the finite size 
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FIG. 1. Schematic view of D and V in the hmit N ^ oo when the synchronization transition at K = Kc is similar to the 
second order phase transition. 



effects in the incoherent regime were addressed in the case of general couphng [ij, |l5| , those in the coherent regime 
are stiU unclear. 

This paper investigates the finite size effects on the long-term fluctuations of the order parameter in the phase 
oscillator model ([1]) by using the diffusion coefficient of the time integral of the order parameter. This diffusion 
coefficient represents the deviation of the time integral of the order parameter from its mean value on the sample 
average. Although this statistical quantity has been used in the large deviation theory, it has not been examined in 
the literature of coupled phase oscillators. Denoting the product of N and the diffusion coefficient as D, we analyze 
the properties of D in the phase oscillator models with the sinusoidal coupling function and also with more general 
coupling functions. 

We show that the scaling law of D with respect to system size N is different for the coherent and incoherent regimes 
in the model ([I} with a general coupling function: D ^ 0{1/N°') with a certain positive constant a in the coherent 
regime, and D ^ 0(1) in the incoherent regime. The scaling law in the coherent regime is anomalous because the 
scaling law of D is ^ 0(1) for both the two regimes in other well-known systems such as the Ising model. Moreover, 
we analytically demonstrate that in the coherent regime, D = in the limit N ^ oo. Therefore, the statistical quantity 
D is useful to qualitatively differentiate between the coherent and incoherent regimes, as illustrated in Fig. [TJ This 
property is not found in other statistical quantities such as the variance of the order parameter. When we denote the 
product of N and this variance by V, it follows V ^ 0{1) with system size N and F 7^ in the limit — >■ cx) in both 
coherent and incoherent regimes Q , as shown in Fig. [Tl 



II. A STATISTICAL QUANTITY D CHARACTERIZING LONG-TERM FLUCTUATIONS OF THE 

ORDER PARAMETER 

We introduce the diffusion coefficient of the time integral of R{t) to characterize the long-term fluctuations of R{t) . 
The variance of the time integral R{s)ds is given as follows: 




t+to 



R{s)ds - {R)tt 



(3) 



where {■)t and [•]s represent the time average over the period from to = to — >■ 00 and the sample average 
of different realizations of ujj which are independently chosen from a certain distribution, respectively. Then, the 
following diffusion law holds: 



D = lim cr^(i)/2t. 



(4) 



which represents the deviation of R{s)ds from its mean value {R)tt on the sample average. 
It should be noted that the statistical quantity D is different from the variance of R{t), given by 



V^N[{{R{to)-{R)t)' 



(5) 



which characterizes instantaneous fluctuations of R{t). For example, if R{t) oscillates periodically, then D = whereas 
V>0. 
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FIG. 2. The time evolutions of the variance <J^{t) of the integrated order parameter (upper) and the correlation function C{t) 
(lower) in the Kuramoto model where h{x) = sin(x) and A'' = 24000. (a)(c) The coherent regime where K = 1.68 > Kc = 
1.59- ■ ■ . (b)(d) The incoherent regime where K = 0.8 < Kc. In both regimes, cr^(t) increases linearly with slope 2Dt after 
a transient period. The correlation function C{t) in the coherent regime has a characteristic time period in which C{t) < 0, 
whereas the form of C{t) in the incoherent regime is almost exponential. Each plot is an average over 30 samples. 



Now, we consider the dependence of D on other statistical quantities. We define the correlation function of R{t) as 
follows: 

C{t) = N[{{R{t + to) - {R)t){R{to) - {R)t))t]s. (6) 

Because C(0) — V, we obtain 

C{t) = Vf{t/r), (7) 

where r and / with /(O) = 1 represent the correlation time and the normalized correlation function, respectively. In 
this paper, r is defined as the value satisfying C(r) = C(0)/2. Then, D satisfies the following equation [l6| : 

1 r°" 

D = - C{s)ds. (8) 



2 „ 

Therefore, from Eqs. (O-®, D can be described using V, r, and / as follows 



D = \vt I f{s)ds. (9) 



III. SCALING LAW OF D WITH SYSTEM SIZE FOR THE KURAMOTO MODEL 

The main result obtained for the scaling law of D is as follows: in the Kuramoto model (Eq. ([T]) with h{x) ~ sin(a;)), 
the asymptotic form of D for large N takes 

^ \0{l/N°-) (coherent regime), 
I 0(1) (incoherent regime), 
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FIG. 3. The scaling property of the diffusion coefficient D with system size in the coherent regime of the Kuramoto model, 
where K = 1.68 > if c = 1.59 ■ • • , and iV = 8000, . . . , 64000. Line fitting yields D ~ iV"^-^^. Each plot is an average over 90 
samples. 
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FIG. 4. The diffusion coefficient D with an increase of system size N in the incoherent regime of the Kuramoto model, where 
K — 0.8 < Kc ~ 1.59 ■ ■ ■ . D fluctuates around a finite value for sufficiently large N. Each plot is an average over 90 samples. 

where a > is a certain constant. 

The cause of the difference in the scaling laws can be intuitively understood from Eq. ([9|) as follows. It is known 
that V ~ 0(1) and r ^ 0{1) for the phase oscillator model ^ and also for other well-known models fT7l.[T8j 

except for the synchronization transition point. If /(s) is written in the simple exponential form, as is commonly the 
case [l3j[l3|) then fis)ds is finite. In fact, the form of /(s) is almost exponential in the incoherent regime of the 
Kuramoto model as numerically shown later. Therefore, with consideration of the above facts, we can infer D ~ 0{1) 
from Eq. ([9]) in the incoherent regime. However, in the coherent regime, /(s) is not in a simple exponential form j^], 
and thereby the scaling law of D is different from D ~ 0{1). 

This section shows the scaling law (|10p numerically. We assume that the distribution of the natural frequencies 
LUj is Gaussian with mean zero and variance one. In numerical simulations, the natural frequencies are generated 
from the Gaussian distribution in a random manner. Figure [5] shows the differences in the time evolutions of cr^(t) in 
Eq. ^ and C{t) in Eq. ([5]) between the coherent and incoherent regimes. The value of D is estimated by fitting the 
values of cr^(t) with a line of slope 2Dt for a sufficiently large t. We separately consider the coherent and incoherent 
regimes. 



A. Coherent regime 

Figure ini shows the dependence of D on the system size N. We clearly see that D is scaled as I? ^ 0{1/N°') where 
a = 1.33 for K — 1.68. In Sec. V, we analytically show that D = in the limit N — oo. 



B. Incoherent regime 

The diffusion coefficient D fiuctuates around a finite value for sufficiently large N as shown in Fig. H) It implies 
D ^ 0(1). We support this scaling property by using Eq. The form of the correlation function C{t) of the order 
parameter is almost exponential as shown in Fig.[SJa) [9]. Therefore, if ^ ~ 0{1) and r ^ 0(1), D should be of 0(1) 
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FIG. 5. The correlation function C{t) for = 64000 (a), the variance V (b), and the correlation time r (c) of the order 
parameter in the incoherent regime of the Kuramoto model, where K = 0.8 < Kc = 1.59---. C{t) almost exponentially 
decreases with t. V and r fluctuate around a flnite value for sufficiently large A'^. Each plot is an average over 90 samples. 

from Eq. (0). In fact, our numerical simulations confirm that V 0(1) and r - 0(1) [1, [H, [H [13, [Hi as 

shown in Figs.ElJb) andEJ^c). Therefore, we conclude D ~ 0(1). 



IV. SCALING LAW OF D WITH SYSTEM SIZE FOR MORE GENERAL COUPLINGS 

This section numerically confirms the scaling law (jlO[) for other representative couplings to enhance the generality 
of our result. The natural frequencies coj are chosen from the Gaussian distribution in a random manner as in the 
previous section. Again we separately consider the coherent and incoherent regimes. 

A. Coherent regime 

In addition to the sinusoidal coupling function treated in the previous section, we consider the following three 
coupling functions [6|: 

(i) A coupling of a generic form with nonzero second harmonic: 
h{x) = sin(2;) - (1/2) sin(2a;), 

(ii) A coupling without the second harmonic term but with the third harmonic term: 
h{x) — sin(a;) — (1/2) sin(3a;), 

(iii) A coupling without the symmetry: 
h{x) — sm{x + 7r/4). 

We choose these coupling functions because most coupling functions are classified on the basis of the value of the 
critical exponent of the order parameter R (for example, see p. 30 and p. 31 in 6]) into the following three cases: the 
coupling function is sinusoidal (Sec. Ill); the coupling function has the second harmonic term ^,3] (case (i)); and the 
coupling function lacks the second harmonic term and possesses the third harmonic term Q (case (ii)). Additionally, 
case (iii) is considered to examine the effect of asymmetry in the coupling function. 

Figures [6ja)-(c) show that D decreases with N as D ^ j-^j, ^^^^^ j^j-j^ ^ jY-i.34 j-^j. ^^^^ ^--j^ ^^^^j jj ^ jY-i.42 

for case (iii). Therefore, the scaling law ([TOl) in the coherent regime holds for these coupling functions. Figures|6l[d)-(f) 
imply D = Q m the limit N oo. In Sec. VI, we analytically show that Z? = in the limit iV — J- oo. 



B. Incoherent regime 

We demonstrate that D ^ 0(1) for any large N in cases (i)-(iii). As shown in Figs. [71^a)-(c), D fluctuates around 
a finite value for sufRciently large N. It implies D ~ 0(1). Further, we support this scaling as follows. The forms 
of the correlation function of the order parameter are almost exponential as shown in Fig s. El^a), El^d), andlHJg). In 
addition, our numerical simulations show that V ~ 0(1) [1, [l3il3 and r ~ 0(1) [TtI. ITq as shown in Figs. [5Jb)(c), 
[5I^e)(f), and [51[h) (i) . Therefore, we conclude D ^ 0(1) from the same argument in Sec. IIIB. 
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FIG. 6. The scaling property of the diffusion coefficient D with system size A'^ (upper) and the correlation function C{t) 
(lower) in the coherent regime of system with (a)(d) h{x) = sin(a;) — (1/2) sin(2a::) and K = 1.75 > Kc = 1.59- ■ • , (b)(e) 
h{x) = sin(3;) - (1/2) sin(3x) and K = 1.675 > Kc = 1.59 ■ • ■ , and (c)(f) h{x) = sm{x + 7r/4) and K = 2.04 > Kc ^ 1.927- • • , 
respectively, where TV = 8000, . . . ,64000 for (a)-(c) and A'^ — 24000 for (d)-(f). Each plot is an average over 90 samples. 
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FIG. 7. The diffusion coefflcient D with an increase of system size N in the incoherent regime of system ([1} with (a) h{x) = 
sin(a; + 7r/4) and K ^ 0.8 < Kc = 1.927- ■■, (b) h{x) = sin(a;) - (1/2) sin(2a;) and K = 1.25 < Kc = 1.59---, and (c) 
h{x) = sin(a;) — (1/2) sin(3a;) and K — 0.8 < Kc = 1.59 • ■ • . D fluctuates around a finite value for sufficiently large A''. Each 
plot is an average over 90 samples. 



V. DERIVATION OF L> = FOR THE KURAMOTO MODEL 



In this section, we analytically demonstrate that D = in the limit iV — )■ oo in the coherent regime of the 
Kuramoto model by reference to Daido which deals with statistical properties in the vicinity of the transition 
point. In general, fluctuations are amplified near the transition point [l3)[l3- Therefore, if I? = is derived in the 
limit K ^ Kc + 0, then D = would be justified for K > Kc- 

We define the complex order parameter Q as follows: 

1 ^ 

where Q represents the frequency of entrainment and \Z{t)\ — R{t). We denote the diffusion coefficient of Z{s)ds 
by Dz, which is defined as follows: 



Dz = hm CT|(t)/2t, 



(12) 
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FIG. 8. The correlation function C{t) for A'^ — 64000 (left), the variance V (middle), and the correlation time r (right) of the 
order parameter in the incoherent regime of system ([!]) with (a)(b)(c) h(x) = sin(a; + 7r/4) and K = 0.8 < Kc ~ 1.927- • • , 
(d)(e)(f) h{x) = sin(a;) - (1/2) sin(2a;) and K = 1.25 < iCc = 1.59---, and (g)(h)(i) h{x) = sin(a;) - (1/2) sin(3a;) and 
K = 0.8 < Kc = 1.59 ■ ■ ■ . C{t) almost exponentially decreases with t. V and r fluctuate around a finite value for sufficiently 
large A^. Each plot is an average over 90 samples. 



where the variance of Jg Z{s)ds is given by 

a|(<) = N 



t+ta 



Z{s)ds - {Z)tt 



to 



With the correlation function of Z{t) given by 

Cz(t) ^ N[{{Z{t + <o) - {Z)t){Z*{h) - {Z*)t))t]s. 
where Z* represents the complex conjugate of Z , Dz satisfies the following equation |16i] : 



Dz = - 
^ 2 



Cz{s)ds. 



In the limit iV — > oo, Cz{t) takes the following form Q: 

Czit) ^ n(yQt)/v/Q, 



(13) 

(14) 
(15) 

(16) 



where n(x) — e^^^^^ +e~^^l'^l/v^~ s|a;|(e~*l^l -j-e"^"!^!) with a certain constant s and Q — Q{K). Integrating the 
above equation by t for t G (— oo, oo), we can obtain Dz = from Eq. (jlSp . 

Next, let us prove that in the limit of — > oo, if Dz = then D = 0. We introduce a variable w{t) for representing 
the fluctuations of Z{t) as follows: 



w(t) = Z{t) 



(17) 



9 



where 

Z = ( hm Z)t, (18) 

A^— >-oo 

and Z = const. ^ in the coherent regime (see p. 29 in Q). For large N, w{t) is small Q in the regime after an initial 
transient period ^2^. Hereafter, let us assume that i = is included in that regime. The deviation w{t) is scaled 
as 0(1/V^) .9, 20]. It should be noted that we can assume Z to be a positive real number because of the rotation 
symmetry of the system [23|. Then, R{t) can be approximated as follows: 

R{t) = ^{Z + Re{w{tW + (Im(u;(t)))2 

w Z^Jl + 2Re{w{t))/Z 

Z + Ke{w{t)), (19) 

where we have neglected the higher-order terms that vanish in the limit N oo. From Eq. (|19p . if the diffusion 
coefficient of Re{w{s))ds is equal to 0, then that of J^, R{s)ds is also equal to because Z is constant. That is, in 
the limit N oo, if Dz — then D — 0. Hence, because Dz = in the limit oo as explained in the previous 

paragraph, _D = in the same limit. 



VI. DERIVATION OF D = FOR A GENERAL COUPLING FUNCTION 



In this section, we analytically demonstrate that Z? = in the limit — > oo in the coherent regime of the phase 
oscillator model ([IJ with a general coupling function h{x) by reference to Daido 9], which considered the sinusoidal 
coupling function. 

Let us explain a key assumption for deriving D = 0. When the system shows synchronization, the oscillators are 
divided into the following two groups: (i) entrained oscillators, which are synchronized with the frequency fi, and (ii) 
nonentrained oscillators, which are not synchronized. Entrained oscillators play a role in reducin g th e fluctuations of 
the order parameter [l3|. On the other hand, nonentrained oscillators show minor fluctuations [ld|- We denote by 
r(w) the distribution of the mean frequencies w of the oscillators. Here, we assume that the following condition holds 
in the coherent regime: 

lim r((Ii) = with ^ oo, (20) 



where il is the common frequency of the entrained oscillators. Daido 21| analytically demonstrated that the above 
condition holds if the coupling function h{x) exhibits only one local minimum and only one local maximum in its 
domain. Further, Daido 5] numerically conflrmed that this condition holds for more general coupling functions. Figure 
ini illustrates a typical example, where Eq. (pn|) holds in a coherent state whereas it does not hold in an incoherent 
state [ij . Equation ([20)1 means that the density of non-entrained oscillators with mean frequencies close to but not 
equal to significantly decreases as w — fi. As a result, fluctuations become minor \ld\ because the more distant the 
mean frequency of an oscillator is from f2, the smaller are the fluctuations caused by the oscillator. 

Considering the power spectrum of VN{R(t) — {R)t), its asymptotic form in the limit — > oo is given by /(w) = 
IZc C{s)e"^'ds dl. From Eq. ©, we obtain 

= i lim /(a.). (21) 

z >-o 

This equation indicates that the value of D is almost determined by the value of r(a;) around uj ~ Cl, since fl can be 
replaced by due to a transformation of the phase variables. Therefore, it is natural that the value of D becomes 
very small if Eq. ([20)) holds. 

This section is organized as follows. First, we transform the original equation ([l} into another expression according 
to Ref. pTj . Next, we derive a self-consistent equation governing the fluctuations of the order parameter. Finally, we 
show Z) = by using the Fourier transform of the self-consistent equation. 
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FIG. 9. Schematic view ofT{Cj). (a) An incoherent state, (b) A coherent state. 



A. Transformation of the phase oscillator model 

The coupling function h in Eq. ^ can be generally represented by the Fourier series as follows: 



(22) 



where qi represents the Zth Fourier coefficient for I — ±1,±2, • • • . We assume that all the synchronized oscillators 
rotate with the common frequency VL (21| | and introduce the generalized complex order parameters as follows: 



1 ^ 



where Zi{t) = Z{t). By using qi and Zj, the order function H{x) [^T*! is defined as follows: 

H{x) = -Y,qiZie-^^''\ 
I 

From Eq. (|24p , we can transform Eq. ([T]) into the following form: 



Oj ^ LUj - KH{ej - nt). 



By introducing new variables 9j — 9j — fit and Aj = tuj — fl, Eq. ()25p is transformed into 

Sj/dt ^ Aj - KHiOj). 



(23) 



(24) 



(25) 



(26) 



B. The self-consistent equation of fluctuations 

First, let us introduce a new variable wi for representing the fiuctuations of Zi(t) as follows: 

wi{t) = Zi{t) - Zu (27) 

where 



Zi = ( lim Zi)t., 



(28) 



and Zi — const. 7^ in the coherent regime (see p. 29 in Q). For large N ^ wi{t) is small [9| in the regime after an 
initial transient period jH]. Hereafter, let us suppose that t = is included in that regime. 
Now we assume that 9j can be divided into two parts as follows @ , 



(29) 
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where 4'j and correspond to the dominant phase motion and the small deviation from it, respectively. The dominant 
phase motion can be described as follows: 

di/jj /dt = Aj - KH{i:j ) , (30) 
V',(O) = 0^j(O) = 0j(O), 

where 

i?(x) = -^g,Zie-2"'^ (31) 
I 

That is, ipj corresponds to 9j in the infinite-size system. 

Next, let us introduce a self-consistent equation of wi{t). If we put 



wi = VNwi, (32) 

then wi is 0(1), as discussed in p. 760 of [9.]. The deviation (pj induced by wi should be of 0{N^^^^), so that we can 
expand 9j in N^'^^^ as follows: 

0;=V^,+</,„ (33) 
= V, + 4= + 0{N-'). 



Substituting Eqs. (H?]), and ([22]) into Eq. (gH) and comparing 0{N^'^/'^) terms, we obtain 

d4>j/dt = if ^ g/(-27ri/Z/0j + wi)e-^'''^'''' . (34) 

Furthermore, from Eqs. (^5]) . ([27]), and we can derive 

N 

wi =VNi-Zi +N-^Y^ e^'^*''^^ ) 

N 

+ 2TrilN-^ J2 ^j-e^"'''''' + 0(A^-^/=^). (35) 

Note that the first term on the right-hand side (r.h.s.) of Eq. ([55]) is 0(1) [o*]. By inserting the solutions of Eq. ([M)) 
into Eq. psp and by considering the limit — >■ oo, we arrive at the self-consistent equations for wi as follows: 

wi{t)^Pi{t)+27TilKy2qi, f dt'Ai,i,{t,t')wi,{t'), (36) 
I' Jo 

where 



TV 

Pi{t)= lini V77(-Z( + 7V-iy e^'^^^'^O, (37) 



and the kernel A;^;/ is defined by 

N 



Au>{t,t')= lim N-^y^exp\2ni{liPj{t)-l'iljj{t')) 

^1 

r'(^,(r))}. (38) 
Here H'{x) represents the derivative of H{x) with respect to x. 



if / drH' 
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C. The Fourier transform of the self-consistent equation 



The goal of this subsection is to show that D — in the hmit — > oo under assumption ((20| . From the 
discussion in Sec. V, Z? = if Dz — 0. Further, the condition Dz = is equivalent to lim^^^Q Iz{i^) = where 
/^(a;) — Cz{s)e^'^'^ds is the asymptotic form of the power spectrum of y/N{Z{t) — {Z)t) in the limit N ^ oo. 
Therefore, it is only necessary to show lim(^_j.o = 0. To evaluate Iz{i^), we cast wi and Pi{t) into the form of 
the frequency domain representation, respectively, as follows (23j : 

wtico) EE hm ^ r dTWi{r)e~"^\ (39) 

and 

Pt[u) EE hm ^ r dTPi{T)e'^--. (40) 
The equation VmY^^^Q Iz{i^) = is satisfied if 



Equation (|4T|) holds if 



limw*(a;)=0. (41) 

a;— 



= lim wUoj) 

= 0, (42) 

which we will show in the rest of this section. 

Let us consider the Fourier transform of both sides of Eq. p6l) : 

w:{Lo) = P,*{uj)+g*{u), (43) 

where g*{'^) represents the Fourier transform of the last term of Eq. ([55]). Concerning the first term of the r.h.s. of 
Eq. (1121), condition ^ yields 

limP;(c.)=0. (44) 

Namely, the first term in the r.h.s. of Eq. ([M)) does not affect the value of D. 

To proceed further, we show that the kernel j4(t,<') is approximately represented in the form of B[xt — yt') with 
certain constants x and y satisfying xy > ii K fa Kc [9]- We divide A as A = + A(ne), where A(o) and 
^(ne) represent the contributions from the entrained and nonentrained oscillators, respectively. Because tpj{t) of the 
entrained oscillators is constant, A(^^j can be represented as follows: 

A(e)(t,t') =S(e)(t-0, (45) 

where ^(o) is a certain function. We can represent A(^^^j as follows: 

A(„e)(i,t') - lim TV-i^e^-^^- (46) 



where 



-^^/T7?'(Vfe(T)). (47) 

In the case of a nonentrained oscillator, if we rewrite Eq. (|T7)) as 

Xk{t, t') ^{l- + Kilt - I't') 

+ Mt,t'). (48) 
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where AJ^, represents the mean frequency of V'fc, then the last term Xk{t, t') is bounded. This is because the last term 
of Eq. (H7)) and ipk{t) are periodic as shown in Appendix A. Because the bounded variation of Xk{t,t') should be 
small compared to the other terms for 1, we can neglect Xk{t,t') as follows [9]: 

Xk{t, t') /')V'fc(O) + A',(/i - I't'). (49) 

This approximation is good enough near the critical point K = Kc [9^] at which the bounded function Xk{t, t') vanishes. 
Therefore, the kernel A^^g) is represented in the following form: 

A(„e)(t,i') =^(nc)("-^'i'), (50) 

where 

= lim iV-iy exp(27ri{(/-/')V'fc(0) + Afe(?i-^Y)})- (51) 

k 

We can exclude the case of < @, in which, for large t and < t* < i, 

= lim N-^Ycxp{2TTi{{l-l')^kiO)+A'J{t+\il'/l)\t*)}) 

k 

= 0. (52) 
From Eqs. (HSI) and ([SO]) , the kernel A can be expressed by using a function B as follows: 

A{t,t') = B{xt-yt'), (53) 

where x{l)y{l') > 0. 

In order to use the Fourier transform, we replace by in the r.h.s. of Eq. ([36]), based on the discussion in 
Appendix B. Then, from Eqs. (pS)) and we obtain 

wUto) = P:{Lo) + 2nilKY,m'^^B*{u:/x)wl{yulx), (54) 

I' 1^1 

where 

/oo 
dTB{T)e-"^\ (55) 
-oo 

Note that we have not divided the r.h.s. of this equation by T. From Eqs. pi)) -([5 ^ . we can obtain 

lim w*i (w) 

= \im2T:ilKy^ qi,^^B*(uj/x)wiAyuj/x). (56) 
tJ-J-O ^-^ \x\ 



Equation ((56)) can be rewritten as follows: 

Ei=Y,Fu,Ev, (57) 

where the coefficient Fi^ii is a certain constant for / = ±1, ±2, • • • and = ±1, ±2, • • • . The trivial solution of Eq. ([57)) 
is E'i = for all I. Let us assume that Eq. ([57)) has another solution Ei = Ei. Then, we can easily show that Ei = cEi 
with any constant c is also a solution of Eq. ([571) . However, this statement contradicts the fact that \Ei\'^ (and D) is 
bounded for K ^ Kc- Therefore, the only solution of Eq. (1551) must be i?; = for all I. 

In fact, in the limit K Kc, we can derive Ei — for all I as follows. Let us divide B*{uj) as B*{uj) = 
B*^^{ijj) + B*^^^-^{u), where B*^^^{lo) and B*^^^-^{uj) represent the contributions from the entrained and nonentrained 
oscillators, respectively. Condition ([20)) yields lim^^o -B^* j^,-) (w) = 0. In the limit K Kc, i?^*^j(aj) = because there 
are no entrained oscillators at isT = Kc- As a result, we obtain Ei = for all because the right-hand sides of 
Eqs. dSD) and ([57]) vanish in the limit K Kc- Consequently, Eq. ([35]) holds. 
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VII. SUMMARY AND DISCUSSION 



We have investigated the statistical properties of long-term fluctuations in the system of globally coupled phase 
oscillators ((T]) with general coupling, by using the statistical quantity D, which is the diffusion coefficient of the 
temporal integration of the order parameter. To understand the finite size effects in the system behavior near 
the synchronization transition point, the scaling property of D with system size N has been examined. We have 
demonstrated that D ~ ©(l/iV") with a certain positive constant a in the coherent regime, and D ~ 0(1) in the 
incoherent regime. The difference in the scaling laws is caused by the difference in the correlations among the phases 
of the oscillators at different times; these correlations remain after a long-term period in the coherent regime. In other 
well-known systems such as the Ising model, the correlation function of an order parameter decays exponentially 
with time [l3, [SIj 8'iid thereby, D follows D ^ 0(1) with respect to the system size N both in the coherent and 
incoherent regimes except for the transition point. For the phase oscillator model ([J), such a difference in the scaling 
laws of D has not been found for other statistical quantities such as the variance and the correlation time of the order 
parameter Q. The scaling property of D in the coherent regime has been further explored in the limit N ^ oo. 
We have analytically demonstrated that D — Q m. the limit iV — >■ cx) for the system with a wide range of general 
coupling functions. If the system exhibits periodic behavior, this result would be trivial. However, this is not the 
case because non-periodic (chaotic) behavior is present even in the coherent regime of the system, as supported by a 
positive Lyapunov exponent [l^]. Although the finite size effects on the statistical properties in the phase oscillator 
model (HI) have been well studied for the sinusoidal coupling function [9l-[l3j , they have remained unclear for a general 
coupling function except for several properties [13, [ill • We have clarified one aspect of the finite size effects in the 
coherent state for coupling functions satisfying Eq. (I^Ul) . which holds for a large class of general coupling functions 

Our result is useful to derive the scaling property of D for the coupling strength interval \K ~ Kc\. From the scaling 
hypothesis [13, [iBl, the variance and the correlation time of the order parameter are scaled a,sV ^ \K — Kc\^'^ and 
T ^ \K — Kc\~^, respectively, where 7 and z are the critical exponents (TTj, JJi]. Combining these scaling laws and 
Eq. ([9]), we can derive the following scaling law: 

/oo 
f{s)ds. (58) 
-00 

From our numerical simulations, we found that, in the limit N — > 00, f{s)ds goes to in the coherent regime 
whereas it is finite in the incoherent regime. Therefore, if the bifurcation of the order parameter is supercritical, we 
obtain the following scaling law in the limit N ^ 00: 

The critical exponent is dependent on 7 and z in the incoherent regime, whereas it is independent of them in the 
coherent regime. It is known that 7 = z = 1 for the sinusoidal coupling function [ol. O [l5l|. 

There are two sources for the order parameter fluctuations. The first is the oscillators that fail to synchronize 
with the order parameter motion. The second is the randomness in the distribution of the natural frequencies. In 
order to show that the main source of the fluctuations is the first one in the coherent regime, we have performed 
numerical simulations by excluding the randomness of the natural frequencies. Namely, the natural frequencies ujj 
are not randomly but deterministically chosen from the Gaussian distribution G(a)) with mean zero and variance one, 
i.e. j/{N + 1) = G{uj)duj. Also in this case, we have obtained the same scaling property of D in the coherent 
regime, i.e. D ^ 0{1/N°-) with a positive constant a for all the coupling schemes considered in this paper. The result 
for the Kuramoto model is shown in Fig. 1101 Confirming the scaling law of D in the incoherent regime should be our 
future work. 
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FIG. 10. The scaling property of the diffusion coefficient D with system size A*' in the coherent regime of the Kuramoto 
model, where K = 1.635 > Kc = 1.59- ■■ , N — 800, . . . ,8000, and ujj is deterministically generated. The line fitting yields 



Each plot is an average over 10 different initial conditions. 



Appendix A 

Because ^pk{t) is periodic, the last term of Eq. (|T7l) is also periodic as follows: 

dTH'(Mr)) 

dr - 

d^k^H'iMr)) 
d^pk 



dipk 



H'i^k 



(Afc - KH{^k)) 



= --\og\{/\k-KH{^k))\ 

= --^log|dVfc(T)/dT|. 



(Al) 



Appendix B 

In order to use the Fourier transform, we consider replacing by in the r.h.s. of Eq. p6|) . First, by defining 
'Wi{t) = for t < 0, we can replace by in the r.h.s. of Eq. For the group of entrained oscillators, 

by defining i?(o)(i) = for t < 0, we can further replace by in the r.h.s. of Eq. (|36|) . For the group of 
nonentrained oscillators, we separately treat the cases oi I = I' and I ^ V . In the case oi I = I' , we can replace J_^^ 
by in the r.h.s. of Eq. ([55)1 by adequately defining = or i?(,ic)(— = for t < for each /. In the case 

of I ^ Z', we replace t' with a large t*{< t) in i3(ne)(i) of Eq. (|5T|) . If (/ — I') ^ Q and both t and t* are sufficiently 
large @, we obtain 



B^ne){lt-l't*) ^ lim N-^"^< 

k 



exp(27rz{(/ - ^OV-fcCO) + Ki^t - I't*)}) = 0. 

iv— >-oo ' — ^ 
k 

As a result, we can replace by in the r.h.s. of Eq. 



(Bl) 
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